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Her  Majesty  the  Queen  in  Right  of  Canada,  as  represented  by  the  Minister  of  National  Defence, 


Part  1  -  Theoretical  Description  of  PCFEM  Code 


PUC:  A  Periodic  Unit  Cell  based  code  for  modeling  materials  with  inclusions 


The  theory  behind  the  implemented  Periodic  Unit  Cell  based  code  is  described.  The  work  is 
based  on  references  [1-3], 
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Figure  1 

The  material  is  periodic  in  the  directions  x  and  y  and  is  excited  by  an  oblique  incidence  plane 
wave  given  by: 


~  /  \  _  j  -y7:(sin#cos^.x+sin0sin^y+cos0z) 

Pi{x,y,z)  —  Ae 


where  we  have  assumed  a  ejat  temporal  dependency. 
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Figure  2 

We  start  with  the  Unit  Cell  Finite  Element  system  (see  Figure  2  for  a  cut  along  plane  xz  and  for 
notations): 


([K]-^[M]  -[c]  V,n 

l  -W 

[u]  and  { p]  contain  the  structural  displacement  (displacements  along  x,y,z)  and  pressure  dof 
respectively  before  imposition  of  boundary  conditions,  {i7}  and  {eft}  represents  the  external 

force  nodal  vector  acting  on  the  structure  and  the  external  normal  pressure  gradient  nodal  vector 
acting  on  the  fluid  domain  boundary  respectively. 

Note  that  for  an  acoustical  excitation  j/7}  =  0  and  {cpj  is  only  non-zero  at  degrees  of  freedom  on 
boundaries  Z+and  2T  .  {<t>}  is  given  by 

{$}=[  {7V(x)}-^-(xWs(x)  (3) 

’  JZ+vjZ~vjBvjTuLuR  <•  V  — ' 

unext 

where  p  is  the  acoustic  pressure  inside  the  finite  element  domain  Q+  u  Q “ .  Actually,  it  will  be 

shown  after  that  only  the  contribution  of  Z+  and  2T  needs  to  be  calculated  (see  Eq(67)). 

Since  the  material  is  periodic  in  direction  x  and  y,  any  function  solution  of  the  problem  must 
satisfy: 

r  (x  +  2dl ,  y  +  ld2 ,  z)  =  r  (x,  v,  z)  e  j2d2qy 

=  r (x,  z) Q1<t>y 

with 


Pr m 


-{O} 


(2) 
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qx=k  sin  0  cos  (j) 
qy  =  k~  smOsm(j) 


(5) 


Note  that  a  function  T  (x,y,z)  satisfying  Eq(4)  can  be  written  as  : 

r(x,y,z)  =  Q(x,y,z)e~jq‘x  e~jqyy  (6) 

where  ®(x  +  2dx,y  +  2d2,z)  =  ®[x,y,z) .  &(x,y,z)  is  a  periodic  function  of  x  and  y  of  period 
2dx  and  2 d2  in  each  direction. 


To  calculate  this  vector,  the  pressure  in  Q£  is  expanded  in  Bloch  series.  The  reflected  pressure 
in  Q;  can  be  written  using  Eq(6)  as: 


Pr{x,y,z) 


&~r(x,y,z) 


e 


-  y'&(sin#cos^x+sin#sin^y) 


(7) 


And  since  0r  (x,y,z)  is  a  periodic  function  in  x  and  y,  it  can  be  expanded  in  term  of  spatial 
Fourier  series  along  x  and  y  as: 


With 


Q;(x,y,z)=  £  ®rmn{z)e,7‘mxe 


m,n=- oo 
+00 


I  ®; 


eiknmz  g-jYi™*  e~jY?ny 


m,n=-co 


0~  = 


_L _ 1_ 

2d,  2 d2 


Jo  'Jo  0,.  (x,  y,  z)  ei7‘mxeJ/'mdxdy 


(8) 


(9) 


and 


Yx 

Yy 


n 

d  i 

n 


d2 


(10) 


For  the  z-dependence,  the  reflected  wave  is  supposed  to  propagate  in  the  -z  direction. 
Therefore  the  total  pressure  in  can  be  written  as: 


p  (x,y,z)  =  pi(x,y,z)+  £  PmneiKnZe 


-  jkmnz  -y'&(sin0cos^x+sin0sin^y)  -jyxmx  -j/yny 


m,n=-cc 

+oo 


(11) 


pi(x,y,z)+  £  pmne]k"mZe-J(Ve 


iP„y 


with 
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am  ~  mYx  +  k  sin  0  cos  (f) 
1 1 3n  =  nyy  +  k  sin  #sin  <j) 

Inserting  Eq(  1 1 )  in  Helmholtz  equation  leads  to 

Similarly,  the  pressure  in  kl+e  can  be  written  as 


p+(x,y,z)  =  £  p+mne 


~ikmn\Z  h+) e~jamX e~jPny 


with 


(*»)2=(*+)2 -*£-/£ 


(12) 


(13) 


(14) 


(15) 


dp  . 


T o  calculate  the  right-hand  side  4T  (see  Eq(67)),  the  idea  is  to  express  —  in  terms  of  the 

1  ’  dz 

nodal  values  of  the  acoustic  pressure  in  the  finite  element  domain. 

The  normal  pressure  gradient  on  Z+  [z.  =  /?+)  can  be  calculated  from  Eq(14)  as 

dp 


dz 


=  V  V  -jk+  p  e  j(ve 

J  mnrmn 


-iPny 


z=h 


m=—M  n=—N 


(16) 


=  (-jK,e-Jve-My){P:,} 

where  ( 2M  + 1)  terms  along  x  and  ( 2N  + 1)  terms  along  y  have  been  kept  in  the  series 
expansion.  The  vector  (-jk*ne~JctmXe~j/3ny^  is  of  size  (2 M  + 1)  x  ( 2N  + 1) . 

The  acoustic  pressure  on  Z+  [z  =  /z+)  can  be  approximated  as 


p+(x,y,h+)=  X  X  "Ve  '' 


jP„y 


(17) 


m=-M  n=—N 


Where  the  coefficient  pmn  writes: 

1  1  c2d<  c2d-. 


P  mn 


2dl  2 d2 


Jo  'Jo  p+{x,y,h+^eictmXeiP"ydxdy 


(18) 


p+(x,y,h+^j  can  be  expressed  in  term  of  pressure  nodal  values  as: 

p+(x,y,h+ )  =  (Vp(x,  y,h+)){p} 


(19) 
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The  vector  ^./V'0  (x,  is  of  size  Np  where  Np  is  the  total  number  of  pressure  dofs. 

Therefore  substituting  Eq(19)  in  Eq(18)  and  putting  the  (2 M  + 1)  x  ( 2N  + 1)  terms  in  a  vector 
leads  to: 


[p~) = xrxrJT  C  (N>  (x.yy^dyW 


2d]  2 d2 


=  [A+]{P} 

Matrix  [.4+J  is  a  rectangular  matrix  of  size  ((2M  +  l)(2./V  +  l),iV'D) 

The  normal  pressure  gradient  on  £“  (z  =  0)  can  be  calculated  from  Eq(1 1 )  as 


dp 


dz 


z= 0 


dp,- 

dz 

dpi 

dz 


M  N 


X  X  jKnPmne~iamXe 


~iP„y 


z=0  m=—M  n=—N 


z= 0 


where  the  same  number  of  terms  as  before  has  been  kept  in  the  expansion. 
The  acoustic  pressure  on  £“(z  =  0)  can  be  approximated  as 


P  (x,  y,  0)  =  p,  (x,  v,  0)  +  X  X  P~mne~i<Ve 


-M,y 


m=—M  n=—N 


(20) 


(21) 


(22) 


Where  the  Fourier  coefficient  p  writes: 

i  mn 


P  mn 


~\0  ‘J0  (p~  (*,  y,  0)  -  A  (x,  y,  ti))ei,veiP"ydxdy 


2d1  2 d2 


(23) 


p  (x,  v,  0)  can  be  expressed  in  term  of  pressure  nodal  values  as: 

^(x,j,0)  =  (a^p(x,j,0)}{/?}  (24) 

pt .  (x,  v,0)  can  also  be  expressed  in  term  of  incident  pressure  evaluated  on  the  FEM  nodes: 

pi{x,y,Q>)  =  (Np(x,y,0)){pl)  (25) 

Therefore  substituting  Eq(24)  in  Eq(23)  and  putting  the  (2 M  + 1)  x  ( 2N  + 1)  terms  in  a  vector 
leads  to: 


'  C  (i V'  (x,y,0))dxdy{p  -  Pi] 

lAd{p-p] 


Finally,  Eq(3)  rewrites  : 


(26) 
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{^}  =  J2+  [N(x,yX)}  { P) 

+L- {#(*>  J'*0)} {-jKne~jamXe  jPny\A  ]ds {p} 

-J2-  {^V(jc,^,0)}  dS  [Pl) 


dS 


dpi 


Besides,  it  can  be  noted  that 

Indeed,  pi(x,y,z)  can  be  expanded  as 

M  N 

P,(*,y,z)=  I 


So  that 


Sp, 

& 


m=—M  n=—N 


M  N 


=  -  y  y  jk  p  eia"xeip"y 

J  mux  i,mn 

Pi.r 


r=0  m=—M  n=—N 

e~iam~e 


=  (-jk;  <,-}** trite 


with 


,  u_L_ Lr"-f”-(e».v« 

i,m«j  0/7  0/7  Jo  J( 


}(A^p(x,  v,0))r/xr/y{y>,} 


2d,  2d2  Jo  Jo 

=  MW 

Therefore  Eq(27)  becomes  : 

' ' °C~ }  =  L  [ N  (*’ y,i h+ )}  {~jk^JamX^jP"y  )[-4+]<7*S'  {/?} 

+|x  {A^(^c,^,0)}  (-jkmneiamXeWny\A-^dS { p } 

dpi 


-2jjN(x,y,0)}fdS 


(27) 


(28) 


(29) 


(30) 


(31) 


(32) 


Or 

(©;j  =  ([A*]  +  [A-])^}-2ji_{w(x,v,°)}^rfS  (33) 

System  (2)  becomes  : 
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(34) 
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Figure  3 


The  nodal  vector  [ u }  and  { p }  can  be  partitioned  in  nine  parts:  surfaces  B,  L,  T,  R,  corner  lines 
Cj ,  C2 ,  C3  and  C4 ,  inner  domain  I  (see  Figure  3). 

{U  }  —  ((wi) (W-ff  )(W-s)(wr)(Mc,  )(mc2  )(mc3  )(mc4  )(w/ ))  (35) 

{p''°,}=((P£)(PJ!>(P2>)(P7-)(Pci)(/»Ci)(^)(Pc4)(P/>)7'  (36) 

Where  (.)  denotes  the  transpose  of  {.} .  {uL),(uR),  (uB),(uT)  contain  the  displacement 
degrees  of  freedom  on  surfaces  L,  R,  B,  T  except  those  on  their  edges.  (uQ  'j ,  (uc^  ^ ,  (uCi  ^ ,  (uCf 
contain  the  displacement  degrees  of  freedom  on  the  unit  cell  edges  along  direction  z.  (uj) 

contains  the  internal  displacement  degrees  of  freedom  (all  those  which  are  not  on  the 
boundaries  of  the  unit  cell).  The  same  notations  are  used  for  the  pressure  degrees  of  freedom. 

Similarly,  the  force  and  normal  pressure  gradient  nodal  vectors  can  be  partitioned  as: 

jF”)  =((Fi>{F,>(FJ)(Fr>(FCi}(FCi)(FCi)(Fc<)(FI))r  (37) 
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{o~"}  =((Oi)<Os><<l>i,>{Or>(®C|){o< 


+ 


O; 


(38) 


K)(o>) 

:®c.){®0(°>) 

where  (<t,i),(Ofi),(Os),((l>7.)  are  nodal  vectors  containing  the  normal  pressure  gradient  to 

surfaces  L,  R,  B,  T1  except  those  on  their  edges.  (|oc  ^Oc  ^Oc  ^Oc  ^  are  corner-line 
nodal  vectors  containing  the  sum  of  the  normal  pressure  gradient  to  surfaces  L  and  B,  L  and  T, 

B  and  R,  L  and  T  respectively2.  (0)  is  the  nodal  vector  which  refers  to  the  PUC  internal  nodes 

not  on  its  lateral  boundaries.  (o>z } ,  (d>z ) ,  (a>z } , (®z ) , (ozQ ),(®c2), (®c3 } > (°c4 ) . (® / )  are  nodal 

vectors  containing  the  normal  pressure  gradient  in  the  z  direction  to  faces  S  and  Z+  .  The  last 
vector  (o)  is  the  nodal  vector  which  refers  to  the  PUC  internal  nodes  not  on  S  and  S+  (i.e 


internal  fluid  nodes  comprised  between  S 

and 

£*)■ 

{®i}=( 

K> 

K)! 

T 

) 

{®«M 

K> 

(®t>! 

)T 

{®»M 
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T 

(39) 


where  the  vectors  have  been  partitioned  in  two  domains  S  ,  S+ . 

=  L,R,B,T,Cl,C2,Ci,C4,I  refer  to  the  vectors  containing  the  normal  pressure 

gradient  in  the  z-direction  on  faces  S  and  S+  . 

There  remains  to  apply  the  boundary  conditions  on  the  unit  cell  which  read  : 


1  Terms  in  — —for  (Oi),(cDfl)  and  in  — ^-for(0B),(0 

'  •  '  '  /-I  t/i  \  /  \ 


dp 


dn 


dn 


2  Terms  in 


dp  dp 
dnx  dn 

a.  y 
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e 


~jVy 


-j\Vx+fy) 


(40) 

(41) 

(42) 

(43) 

(44) 


The  same  holds  for  the  pressure  dots.  These  equations  allows  for  expressing  the  nodal  vectors 
in  terms  of  reduced  nodal  vectors  |wcj  =  and 

{pc}  =  {{pl}{pb){pq){pi))T  ■ 


With 


and 


{«"■}= [e"]{«i 

{p™}  =  [Q"]{p‘} 


[e"]= 


[Q“]  [°1  [°] 

[o]  [e-]  [o] 

[0]  [0]  [e* 

[0]  [0]  [»] 


[0] 

[0] 

[0] 

M 


(45) 

(46) 


(47) 


r 

1c , 


-J\9,+<Py) 


(48) 

(49) 


(50) 
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In  the  previous  equations  [/"]  denotes  the  identity  matrix  of  size  equal  to  the  number  of 
structural  degrees  of  freedom  belonging  to  entity  A  (face  or  edge).  Similarly 


[fi']- 


[ea]  [o]  [o] 

[o]  [2'-]  [o] 

[o]  [o]  [e* 

[0]  [0]  [0] 


[0] 

[0] 

[0] 

[>'] 


(51) 


and 


(52) 


(53) 


(54) 


For  the  external  forces,  we  have  the  following  relationships: 


{^}  +  {^}^={0} 

(55) 

{Fr}  +  {F„je-*={0} 

(56) 

{FC| }  +  {Fc,  j  e»-  +  {Fc>  j  e»’  +  {Fc,  j  <f  ^  =  {0} 

(57) 

Let  assume  that  each  component  of  corner  line  nodal  vectors  |fc  j ,  j 

proportional  to  the  corresponding  component  of  |fc  j ,  that  is: 

and  {FQ}  is 

FC2,i  =  ™iFCui’  FC„i  =  $iFCui’  FCi,i  =  ZiFCui 

(58) 

Then  using  Eq(57),  we  get  the  following  relationship  between  the  proportionality  constants 

<9,  and 


(59) 


1  +  rofe"M  +  &fejVy  +  =  {O} 

The  same  relations  Eq(55)  to  Eq  (59)  holds  for  the  normal  pressure  gradients  but  with  different 
constants. 


We  then  have 


and 


[e1= 


Fren]  =  [_Qf]{fc 

'[e*]  [o]  [o] 

[o]  [e'-]  [o] 

[o]  [o]  [e* 


[0] 

[0] 

[0] 


[o]  [o]  [o]  [/;] 


[efi]= 


[s' 


■M* 

-['»] 


e 


lC, 


X+(Pv 


(60) 


(61) 


(62) 


(63) 


(64) 


_[&yw 

Matrices  [<9q]  anc*  [£q]  are  diagonal  matrices  of  size  the  number  of  structural  dots  on 

corner  line  C,  whose  generic  term  are  er''.  ..  =  arf  ,  ..  =  5/'  „  =  4f  respectively. 

We  also  have 

{<t”}=[e*]{®'’}+{®.}  (65) 

with 

{«*'} =((®i)(®»)K  )(®-»r  <66> 

and 
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(67) 


o>  =  on  o 


07  01  O;  o 


°;>(°> 


=((®:>(o> 


For  the  normal  gradients  in  the  plane  (x,y)  we  have, 


and 


[e*]= 


[e“]  [o] 

[o]  [e*-] 


[o] 

[0] 


[o] 

[0] 


[«]  [0]  ' 
[0]  [0] 
[e-]  [o] 
[0]  [/,«] 


o 

— 

;< 

e1*’ 

■jyPx+Vy ) 


(68) 


(69) 


(70) 


(71) 


Matrices  ,  [i9®]  and  J  are  diagonal  matrices  of  size  the  number  of  pressure  dofs  on 
cornerline  C,  whose  generic  term  are  gj®  ..  =  gj®  ,  5®  u  =  3f  £®  „  =  ^®  respectively. 

Matrices  [/® ] ,  [/®  ] ,  [y®  J  are  diagonal  matrices  of  size  the  number  of  pressure  dofs  on  the 
respective  partition. 

Substituting  Eq(45),  Eq(46),  Eq(60)  and  Eq(65)  in  Eq(2),  multiplying  the  first  and  the  second  line 
of  the  resulting  system  by  [Qu~\  and  [<2P]  respectively,  the  following  reduced  system  is 
obtained: 
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f  [Z“]  [Z"PJ 

{{«■} 

'I 

[[z-]‘r  [Z-] J 

IM 

^?([eT[e*]Kh[e'r  {*■})' 

[z"]  =  [2"]‘r[[K]-®J[M]][e“] 


[z'HC'f 


[z-]  =  -[2“f[C][2'] 


Using  Eq(47)  and  Eq(61)  on  the  one  hand  and  Eq  (51)  and  Eq(68)  on  the  other  hand,  it  can  be 
shown  that: 


[e-JVJM 


S  {0} 

no  internal  force 
acting  on  the  structure 


(73) 


and 


[g'TV]{®‘ 


{0} 

{0} 

{0} 

11°}. 


(74) 


Eq(72)  becomes 


( [zm]  \zup^[ 

M 

T 

'{0} 

r^1 

* 

rN_l 

§ 

< 

M 

>  = 

p,co  L  J  1  ’ 

^  J  > 

(75) 


System  in  Eq(75)  is  Hermitian  and  can  be  solved  using  appropriate  resolution  algorithm.  The 

l"l 


matrices  [IK]  -  or  [M]] ,  [c]and 


CO" 


■[e] 


are  obtained  directly  from  Novafem  without 


imposing  any  boundary  conditions  on  the  fluid  and  the  structure.  System  (75)  requires  to  build 
the  condensation  matrices  [_Q“~\  and  [£?P]-  Note  that  in  jO;|  only  the  normal  incident  pressure 

gradient  degrees  of  freedom  (on  2T  (z  =  0) )  are  non-zero. 
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The  steps  in  the  methodology  are  the  following: 

1 .  Identify  corner  line,  lateral  face  and  internal  nodes 

2.  Partition  global  solution  vector  into  nine  components 

3.  Build  the  periodic  boundary  condition  condensation  matrix  for  the  structure  [ _Q“~\ 

4.  Build  the  periodic  boundary  condition  condensation  matrices  for  the  fluid  \_QP~\  and 

[el 

5.  Calculate  matrices  |^A+J  and  [a~] 

6.  Calculate  the  normal  incident  pressure  gradient  nodal  vector  [O] 

7.  Build  system  (72)  (Hermitian  projection  of  Eq(2)  on  \j) 

8.  Solve  the  projected  Hermitian  system 


1- 

[2“] 

[°] 

J“ 

[°] 

[e'l 

Ce-Jve~J^ 


Calculation  of  |^A+J 
|^A+J  can  be  rewritten  as 

where  [Z)+J  is  a  diagonal  matrix  of  dimension  ((2M  + 1)  x  (2N  + 1) , (2 M  + 1)  x  (2N  + 1)) 

0 


(76) 


-Km-n  0  0  0 

o  -jklM-N+ 1  0  0 

0  0  0 

0  00  -jk 

0  0  0  0 


MN-l 


0 

0 

0 

0 

-A 


MN . 


(77) 


where  [^4+1]  is  a  matrix  of  dimension  y{2M  +  l)  x  (27V +  l), Npj  given  by 
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(78) 


[^+1]  =  Jr  [e^e^]  (Np  (x,y,h+))dS 

=  [eJ'a^eMy‘  ]js+  {iVp  (x,  v,/?+)}  (Np  [x,y,h+))dS 

=  [e;<vto][cY] 

[Y]  =  is  a  matrix  of  dimension  ((2M  +  l)  x  (2iV  + 1), Np)  which  contains  on  each 

line  the  nodal  values  of  term  eja,"xejPny  .  Actually,  only  the  dof  on  face  Z+  are  considered  in  the 
calculation. 


Calculation  of  the  transmitted  power  (power  exchanged  between  the  structure  and  Q+) 

(79) 


W+ =-y\  |~f  p+v.ndS 
2  [_Jan+^ 


Calculation  of  the  power  exchanged  between  the  structure  and  Q 

W-  =  -m\{  pv.ndS 
2  LJar^  J 

=  ^[-j(o(p-)[C]{u} 


(80) 


Calculation  of  the  Reflection  and  transmission  coefficients 


In  the  case  where  only  the  normal  modes  taken  care  of,  the  reflection  coefficient  is  given  by: 


1  1  r2di  c2d2  , 

R  =  P00=— — J0  J0  P  (YY, 0)e 


2dl  2 d2 


jk  sin# cos <f)x ^jk  sin  6  sin  (f> 


’  dxdy 


1  1 


jk  sin  6  cos  jfXj  jk  sin  0  sin  (j>yi 


[C']M 


(81) 


2d1  2  d2 

Similarly,  the  transmission  coefficient  is  given  by  p^0  with: 
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t=k= V  p'  (*.  *  ■ h' ) 


2dt  2d 2 

1  1 


2dx  2 d2 


jk  sin  6 cos  t/>Xj  jk  sin  0 sin  <j>yi 


(82) 


To  validate  these  statements,  note  that  the  incident  pressure  pf(x,y,z)  can  be  expanded  as 

M  N 

Pi(x,y,z)=  X  X  Pumne  jk™ze  ja"xew"y  (83) 

with 


m=—M  n=—N 


I  P,„ )  =  jfC  C  { e‘°-" )  (JV'  (x,  y,  0))  dxdy  { p, } 


2d1  2 d2  Jo  Jo 

=M  {Pi} 

Compare  to  the  reflected  pressure: 


(84) 


p(x,y,z)=  X  pmne'k^e  ja-xe 


iP„y 


m.n— — co 


At  z=0: 


M  N 


Pi(x,y,  0)=  X  X  Pi,mne~Jan,Xe~M> 


m=—M  n=—N 


and  thus. 


P  (x,  y,  0)  =  X  Plne  iamXe  ,P"y  =  RPi  (x,  y,  0) 


m.n—— oo 


(85) 


(86) 


(87) 


Explicitly, 

1  1  f2rf,  f2</;  +°° 


2dx  2 d2 


W,/7=— 00 


1  1  f2</,  f2</: 


(88) 


2dl  2 d2 


n*  i 


m=—M  n=—N 


So  that, 
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(89) 


Pin  =  RPi,,nn  =>  R  =  ~  *  /0>  »)  =  =  P* O' 


Pu 


Pi. 


00 


Note  that  we  can  verify  : 


J _ 1_ 

2dl  2 d2 


Prnn  = 


1  1  r2di  r2rf: 


2dl  2 d2 


Jo  ‘Jo  ^  ( x’  -v’ 0 )  ( )  dxdy  : 


R——r\  2diA 

2d !  2J,  Jo  Jo 


PK  an  \ 
j~rx  i-ry 

a.  d1 

e  1  e  2 


v 


dxdy  =  /?/f  =  /?/),.  ( 


y 


The  same  proof  can  be  given  for  the  transmission  coefficient. 
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Part  2  -  Validation  of  PCFEM  Code: 
Materials  with  no  Inclusions 
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Part  3  -  Validation  of  PCFEM  Code: 
Materials  with  Inclusions 
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PCFEM  -  Validation  cases  with  inclusions 


CASE  1  :  DOUBLY  PERIODIC  CICURLAR  CYLINDRICAL  AIR  INCLUSIONS 


REFERENCES 

•  Easwaran,  V.  and  Munjal,  M.  (1993),  Analysis  of  reflection  characteristics  of  a  normal  incidence  plane  wave  on 
resonant  sound  absorbers:  A  finite  element  approach,  The  Journal  of  the  Acoustical  Society  of  America,  93, 

1308. 

•  Hladky-Hennion,  A.C.  and  Decarpigny,  J.N.  (1991),  Analysis  of  the  scattering  of  a  plane  acoustic  wave  by  a  doubly 
periodic  structure  using  the  finite  element  method:  Application  to  Alberich  anechoic  coatings,  The  Journal  of  the 
Acoustical  Society  of  America,  90,  3356. 

DESCRIPTION 

•  Alberich  anechoic  coating  with  circular  cylindrical  air  inclusion  immersed  in  water  on  both  sides. 

•  Normal  incidence  acoustic  wave. 

•  Air  is  not  modelled  in  References  -  see  as  void. 

•  Only  transmission  results  are  available. 

PCF  MODEL 

•  \PCFEM\work\Test  cases  with  inclusion\casel\casel.pcf 

•  Waveguide  solver  with  TETRA  10  elements 

•  V4  criterion  CPU  time  is  19  s;  Number  of  elements:  1  268  Number  of  nodes:  2  474 

•  V6  criterion  CPU  time  is  98  s;  Number  of  elements:  5  080  Number  of  nodes:  8  784 


Elastic  properties  of  coating  material 

E  (Pa) 

1.4x10s 

v(-) 

0.49 

%  (%) 

23 

p  (kg/m3) 

1100 
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PCFEM  -  Validation  cases  with  inclusions 


W  Figure  1:  Domain  Mesh 
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□  Figure  2:  Transmission  Coefficient  for  casel.pcf 
File  Edit  View  Insert  Tools  Desktop  Window  Help 


r^l  a  IteaJ 


d  a  h  &  |  fc  \%  %  ©  ®  sffi  /.  -  |  a  |  □  m  ■  o 


Good  correlation  obtained  compared  with  reference. 


CASE  2  :  DOUBLY  PERIODIC  CICURLAR  CYLINDRICAL  AIR  INCLUSIONS 


REFERENCES 

•  Easwaran,  V.  and  Munjal,  M.  (1993),  Analysis  of  reflection  characteristics  of  a  normal  incidence  plane  wave  on 
resonant  sound  absorbers:  A  finite  element  approach,  The  Journal  of  the  Acoustical  Society  of  America,  93, 

1308. 

•  Hladky-Hennion,  A.C.  and  Decarpigny,  J.N.  (1991),  Analysis  of  the  scattering  of  a  plane  acoustic  wave  by  a  doubly 
periodic  structure  using  the  finite  element  method:  Application  to  Alberich  anechoic  coatings,  The  Journal  of  the 
Acoustical  Society  of  America,  90,  3356. 

DESCRIPTION 

•  Alberich  anechoic  silicon  coating  with  circular  cylindrical  air  inclusion  immersed  in  water  on  both  sides. 

•  Normal  incidence  acoustic  wave. 

•  Air  is  not  modelled  in  References  -  see  as  void. 

PCF  MODEL 

•  \PCFEM\work\Test  cases  with  inclusion\case2\case2.pcf 

•  Waveguide  solver  with  TETRA  10  elements 

•  'k/l  criterion  CPU  time  is  95  s  Number  of  elements:  6  314  Number  of  nodes:  9  888 

•  X/3  criterion  CPU  time  is  581  s  Number  of  elements:  21  674  Number  of  nodes:  32  240 

•  V3  yields  results  very  very  close  to  X/2  for  a  fraction  of  time. 
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PCFEM  -  Validation  cases  with  inclusions 
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Good  correlation  obtained  compared  with  FEM  and  EXP  results  from  reference. 
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Validation  cases  with  inclusions 


CASE  3  :  SINGLE  PERIODIC  INFINITE  RECTANGULAR  CYLINDRICAL  AIR  INCLUSIONS 


REFERENCES 

•  Hladky-Hennion,  A.C.  and  Decarpigny,  J.N.  (1991),  Analysis  of  the  scattering  of  a  plane  acoustic  wave  by  a  doubly 
periodic  structure  using  the  finite  element  method:  Application  to  Alberich  anechoic  coatings,  The  Journal  of  the 
Acoustical  Society  of  America,  90,  3356. 

DESCRIPTION 

•  Alberich  anechoic  polyurethane  coating  with  an  infinite  rectangular  air  inclusion  immersed  in  water  on  both 
sides. 

•  It  is  similar  to  a  2D  case  (arbitrary  1  cm  height  is  imposed). 

•  Normal  incidence  acoustic  wave. 

•  Air  is  not  modelled  in  References  -  see  as  void. 


PCF  MODEL 

•  \PCFEM\work\Test  cases  with  inclusion\case3\case3.pcf 

•  Waveguide  solver  with  TETRA  10  elements 

•  V4  criterion  CPU  time  is  12  s  Number  of  elements:  720  Number  of  nodes:  1  440 

•  V6  criterion  CPU  time  is  38  s  Number  of  elements:  2  412  Number  of  nodes:  4  228 


Elastic  properties  of  polyurethane  material 

E  (Pa) 

2.81x10s 

v  (-) 

0.479 

m  (%) 

45 

4s  (%) 

1.78 

p  (kg/m3) 

1100  (estimation) 

Mecanum  inc. 
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PCFEM  -  Validation  cases  with  inclusions 


H  Figure  1:  Transmission  Coefficient  for  case3.pcf 


,  |  B 


Good  correlation  obtained  compared  with  reference  except  above  15  000  Hz.  We  have  checked  several  meshes,  and  the 
curve  doesn't  change.  However  the  results  are  very  sensitive  to  damping.  Here  we  apply,  separate  damping  to  E  and 
Poisson  and  deduce  G...  So  we  are  wondering  if  there  is  an  error  in  the  data  especially  the  damping... It  is  surprising  that 

only  case  1  was  reproduced  in  various  papers. 
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Validation  cases  with  inclusions 


CASE  4  :  DOUBLY  PERIODIC  ALBERICH  ANECHOIC  COATING  ON  STEEL  PLATE  -  WATER/SLAB/AIR 


REFERENCES 


•  Meng,  H.,  Wen,  J.,  Zhao,  H.,  Lv,  L.  and  Wen,  X.  (2012),  Analysis  of  absorption  performances  of  anechoic  layers 
with  steel  plate  backing,  The  Journal  of  the  Acoustical  Society  of  America,  132,  69. 

DESCRIPTION 

•  Alberich  anechoic  coating  with  a  cylindrical  air  inclusion  as  in  CASE  1  backed  by  a  STEEL  plate. 

•  Water  on  incident  side  and  air  behind  the  steel  plate. 

•  Normal  incidence  acoustic  wave. 


PCF  MODEL 

•  \PCFEM\work\Test  cases  with  inclusion\case4\case4.pcf 

•  Waveguide  solver  with  TETRA  10  elements 

•  V3  criterion  CPU  time  is  47  s  Number  of  elements:  1  050  Number  of  nodes:  1  996 

•  V6  criterion  CPU  time  is  539  s  Number  of  elements:  9  452  Number  of  nodes:  15  310 


Elastic  properties  of  material 

Rubber 

Steel 

Water 

E  (Pa) 

1.4e8 

2.1ell 

n/a 

v  (-) 

0.49 

0.3 

n/a 

He  (%) 

0.23 

n/a 

n/a 

p  (kg/m3) 

1100 

7890 

1000 

c  (m/s) 

1489 
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|  Figure  7:  Domain  Mesh 
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a  Figure  11:  Absorption  Coefficient  for  case4.pcf 
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PCFEM  -  Validation  cases  with  inclusions 


Good  correlation  with  reference.  However  we  have  tried  to  explain  differences.  We  have  checked  several  meshes  and 
also  meshed  various  air-layers  upstream  and  downstream  slab...  the  curves  do  not  change.  We  believe  that  our  results 
are  more  representative  compared  the  reference  since  (1)  the  anechoic  coating  alone  is  exactly  case  1  that  was 
validated;  and  (2)  the  homogeneous  case  (=without  the  inclusion)  agrees  with  the  Transfer  Matrix  Method  (TMM)  -  see 

next  page. 

Still,  we  are  trying  to  find  out  what's  wrong  with  the  input  data! 


Homogeneous  case  (without  inclusion) 

Comparison  between  PCFEM  (waveguide  solver  * )  and  TMM  (green  curve) 


Mecanum  inc. 


RP  2014-03-30  21:25 
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PCFEM  -  Validation  cases  with  inclusions 


CASE  4B  :  DOUBLY  PERIODIC  ALBERICH  ANECHOIC  COATING  ON  WATER  -  WATER/SLAB/WATER 


REFERENCES 

•  Meng,  H.,  Wen,  J.,  Zhao,  H.,  Lv,  L.  and  Wen,  X.  (2012),  Analysis  of  absorption  performances  of  anechoic  layers 
with  steel  plate  backing,  The  Journal  of  the  Acoustical  Society  of  America,  132,  69. 

DESCRIPTION 

•  Alberich  anechoic  coating  with  a  cylindrical  air  inclusion  with  water  of  both  sides  =  CASE  1. 

•  As  in  case  1,  however  this  time  in  absorption. 

•  Normal  incidence  acoustic  wave. 


PCF  MODEL 

•  \PCFEM\work\Test  cases  with  inclusion\case4b\case4b.pcf 

•  Waveguide  solver  with  TETRA  10  elements 

•  A./6  criterion  CPU  time  is  176  s  Number  of  elements:  6  410  Number  of  nodes:  10  936 
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H  :igure  8:  Domain  Mesh 
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□  Figure  3:  Dissipation  Coefficient  for  case4b.pcf 


is 


File  Edit  View  Insert  Tools  Desktop  Window  Help 


k  i  \ \ o ® « d -  a i  □  m 


PCFEM  -  WG  (a) 
PCFEM -WG  (ad|ss) 
■  Reference 


4000  6000 

Frequency  (Hz) 


Mecanum  inc. 


RP  2014-03-30  21:25 


39 


PCFEM  -  Validation  cases  with  inclusions 

Good  correlation  with  reference  with  same  conclusions  as  in  case  4  relatively  to  the  differences. 

—  I  |2 

NOTE  that  usually  the  absorption  coefficient  is  defined  as  a  =  1  - \R\  .  Here  the  absorbance  they  defined  seems  to  be 

equal  to  adiss  —  D  =  1  -  \R\  -  7”  which  is  in  fact  the  energy  dissipated  in  the  coating.  While  for  the  previous  case 
both  absorptions  are  very  close  due  to  the  steel  plate  causing  very  low  value  of  T,  in  this  case  this  is  no  longer  true. 


CASE  5  :  DOUBLY  PERIODIC  COATED  RIGID  SPHERE  COATING  ON  STEEL  PLATE  OR  ON  WATER 


REFERENCES 

•  Meng,  H.,  Wen,  J.,  Zhao,  H.,  Lv,  L.  and  Wen,  X.  (2012),  Analysis  of  absorption  performances  of  anechoic  layers 
with  steel  plate  backing,  The  Journal  of  the  Acoustical  Society  of  America,  132,  69.  (FIGURE  6) 

DESCRIPTION 

•  Doubly  periodic  ALUMINUM  coated  sphere  in  a  host  rubber  backed  by  a  STEEL  plate  or  WATER.  Water  on 
incident  side  and  air  behind  the  steel  plate  if  not  water  backing.  The  core  of  the  coated  sphere  is  in  Aluminum, 
while  is  coating  is  a  soft  silicon  rubber. 

•  Normal  incidence  acoustic  wave. 


PCF  MODEL 

•  \PCFEM\work\Test  cases  with  inclusion\case5\case5.pcf  (for  steel  backing) 

•  \PCFEM\work\Test  cases  with  inclusion\case5b\case5b.pcf  (for  steel  backing) 

•  Waveguide  solver  with  TETRA  10  elements 

•  A./10  criterion  CPU  time  is  2202  s  Number  of  elements:  29  184  (24  576)  Number  of  nodes:  42  471  (35 
937) 

•  X/6  criterion  CPU  time  is  115  s  Number  of  elements:  6  528  (3  456)  Number  of  nodes:  10  115  (5  491) 


Mecanum  inc. 
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PCFEM  -  Validation  cases  with  inclusions 


VlO  criterion  for  layer  1.  Manual  for  layer  2. 


□  Figure  5:  Dissipation  Coefficient  for  case5.pcf  r^irwis&ii 
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Good  correlation  with  reference.  Finer  mesh  yields  better  correlation  around  peaks. 
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PCFEM  -  Validation  cases  with  inclusions 


CASE  6  :  DOUBLY  PERIODIC  COATED  SPHERE  COATING  -  WATER/SLAB/AIR 


REFERENCES 


•  Wen,  J.,  Zhao,  H.,  Lv,  L.,  Yuan,  B.,  Wang,  G.  and  Wen,  X.  (2011),  Effects  of  locally  resonant  modes  on  underwater 
sound  absorption  in  viscoelastic  materials,  The  Journal  of  the  Acoustical  Society  of  America,  130,  1201.  (FIGURE 

31 

DESCRIPTION 

•  Doubly  periodic  STEEL  coated  sphere  in  a  host  rubber  backed  by  a  STEEL  plate.  Water  on  incident  side  and  air 
behind  the  steel  plate.  The  core  of  the  coated  sphere  is  in  Aluminum,  while  is  coating  is  a  soft  silicon  rubber. 

•  Normal  incidence  acoustic  wave. 


PCF  MODEL 

•  \PCFEM\work\Test  cases  with  inclusion\case6\case6.pcf  (for  steel  backing) 

•  Waveguide  solver  with  TETRA  10  elements 

•  V8  criterion  CPU  time  is  901  s  Number  of  elements:  18  816  Number  of  nodes:  27  753 

•  V10  criterion  CPU  time  is  s  Number  of  elements:  34  992  Number  of  nodes:  50  653 
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□  Figure  6:  Dissipation  Coefficient  for  case6,pcf 
File  Edit  View  Insert  Tools  Desktop  Window  Help 


I  °  II  B 


Frequency  (Hz) 


Mecanum  inc. 


RP  2014-03-30  21:25 


42 


PCFEM  -  Validation  cases  with  inclusions 

Good  correlation  with  reference;  however  some  discrepancies  are  not  explained.  We  thought  that  coarser  mesh  was 
not  enough  to  discretize  geometry;  however  coarse  and  very  fine  meshes  yield  same  results.  In  the  paper,  only  sound 
speeds  in  media  were  given.  So  we  deduced  from  equations  the  properties  given  in  the  table.  We  expect  the 
discrepancies  result  from  errors  in  the  input  data  given  in  the  reference  paper. 
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CASE  7  :  RUBBER/STEEL  PLATE  WITH  DOUBLY  PERIODIC  PRISMATIC  VOIDS  AT  JUNCTION  (PUC  CASE) 


REFERENCES 

•  Email  by  Jeff  Szabo  entitled  "Test  Case"  and  sent  on  February  20th ,  2014. 

DESCRIPTION 

•  Rubber/Steel  plate  with  doubly  periodic  prismatic  voids  at  junction. 

•  Normal  incidence  acoustic  wave. 

•  Properties  of  material  are  frequency  dependent  and  given  in  MAT  file  "rubber3.mat"  and  "steel. mat"  of 
reference. 

PCF  MODEL 


•  \PCFEM\work\Test  cases  with  inclusion\case7\case7.pcf  (for  lxl  PUC) 

•  \PCFEM\work\Test  cases  with  inclusion\case7a\case7a.pcf  (for  4x4  PUC) 

•  Waveguide  solver  with  TETRA10  elements 

•  ~\/2  (USER  DEFINED)  criterion  CPU  time  is  60  s  (347  s)  Num.  of  elements:  1  008  (16  128) 

Num.  of  nodes:  1  836  (24  615) 

•  Dynamic  complex  material  properties.  Their  spectra  are  stored  in  MAT  files:  Rubber3_table.mat  and 
Steel_table.mat  located  in  folder:  \PCFEM\work\Material  Library 


x 


A  y 

1.5  cm  2.13  2.54 


hx=hy=[2.5  5  2.5]  and  hz=[2.5  2.5  8] 
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PCFEM  -  Validation  cases  with  inclusions 


□  Figure  2:  Domain  Mesh 


hx=hy=  [2.5  5  2.5  5  2.5  5  2.5  5  2.5]  and  hz=[2.5  2.5 

8] 


For  the  same  mesh  definition,  one  lxl  PUC  cell  yields  exactly  the  same  results  as  one  4x4  PUC  cell 
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